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Abstract 

One of the main approaches for modeling fracture and crack prop- 
agation in sohd materials is adaptive insertion of cohesive elements, in 
which line-like (2D) or surface-like (3D) elements are inserted into the 
finite element mesh to model the nucleation and propagation of failure 
surfaces. In this approach, however, cracks are forced to propagate along 
element boundaries, following paths that in general require more energy 
per unit crack extension (greater driving forces) than those followed in 
the original continuum, which in turn leads to erroneous solutions. In 
this work we illustrate how the introduction of a discretization produces 
two undesired effects, which we term mesh-induced anisotropy and mesh- 
induced toughness. Subsequently, we analyze those effects through polar 
plots of the path deviation ratio (a measure of the ability of a mesh to 
represent straight lines) for commonly adopted meshes. Finally, we pro- 
pose to reduce those effects through K-means meshes and through a new 
type of mesh, which we term conjugate-directions mesh. The behavior 
of all meshes under consideration as the mesh size is reduced is analyzed 
through a numerical study of convergence. 

1 Introduction 

The classical cohesive theory of fracture finds its origins in the pioneering works 
by Dugdale, Barenblatt and Rice [HOIS]. In their work, fracture is regarded as 
a progressive phenomenon in which crack formation takes place across a cohe- 
sive zone ahead of the crack tip and is resisted by cohesive tractions. Cohesive 
zone models are widely adopted by scientists and engineers perhaps due to their 
straightforward implementation within the traditional finite element formula- 
tion. Some of the mainstream technologies proposed to introduce the cohesive 
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theory of fracture into finite element analysis are the eXtended Finite Element 
Method (X-FEM) and cohesive elements. 

Sukumar et al. [3] first utilized the X-FEM for modeling 3D crack growth by 
adding a discontinuous function and the asymptotic crack tip field to the finite 
elements. Subsequently, the method was extended to account for cohesive cracks 
[S]. It is worth noting that while the X-FEM approach can potentially deal 
with arbitrary crack paths, it becomes increasingly complicated for problems 
involving pervasive fracture and fragmentation. 

On the other hand, the cohesive element approach consists on the insertion 
of cohesive finite elements along the edges or faces of the 2D or 3D mesh cor- 
respondingly [3 [71 [HI IS] . Even though this approach is well suited for problems 
involving pre-defined crack directions, a number of known issues affect the its 
accuracy when dealing with simulations including arbitrary crack paths, namely, 

(i) problems with the propagation of elastic stress waves (artificial compliance), 

(ii) spurious crack tip speed effects (lift-off), and (iii) mesh dependent effects 
(c.f. [10) for a comprehensive review). Despite these well known limitations, 
the robustness of the method makes it one of the most common approaches for 
pervasive fracture and fragmentation analysis. 

Some of the limitations present in early approaches to cohesive element mod- 
els were addressed by subsequent research efforts. For example, artificial com- 
pliance and lift-off effects can be avoided by using an initially rigid cohesive 
law [5] or, more elegantly, a discontinuous Galerkin formulation with an acti- 
vation criterion for cohesive elements jlll I12| . However, the problem of mesh 
dependency is still an active area of research. 

Mesh-dependent effects are direct consequence of cracks being able to propa- 
gate only across boundaries between bulk finite elements. That is, the topology 
of the mesh forces cracks to follow paths that in general require more energy per 
unit crack extension (greater driving forces) than those followed in the original 
continuum. In this work, we first focus on the effects that common topologies 
have on two main mesh-dependent effects, namely: mesh-induced toughness 
and mesh-induced anisotropy. We then illustrate how to decrease mesh-induced 
anisotropy through K-means meshes. Finally, we introduce a new type of mesh, 
termed conjugate-directions mesh, which greatly alleviate both effects. 

The reminder of the paper is organized as follows: in section [2] we intro- 
duce a few preliminary concepts and discuss the mesh dependent effects; in sec- 
tion [3l we define K-means and conjugate-directions meshes and we study their 
mesh-dependent behavior, including a numerical convergence analysis; finally, 
in section [31 we summarize the main conclusions of the article. 

2 Mesh- induced effects 
2.1 Basic definitions 

Let us consider a 2-dimensional brittle fracture problem in which an arbitrary 
crack of length Lc develops, as shown in Figure [1^. The fracture energy Ef 
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Figure 1: (a) Curve representing an arbitrary crack (solid line) and the corresponding 
rectification (dashed line); (b) shortest euclidean distance and shortest path over the 
mesh between two random points. 

required to generate such a crack is proportional to the crack length Lc, that is, 



where Gc is the energy release rate of the material. This implies that, even 
assuming a correct model for the energy release rate of the material, the ac- 
curacy on any prediction of the fracture energy directly depends on a correct 
computation of the crack path length. 

In the most general case, the crack path length could be computed by rec- 
tification of the irregular curve, i.e., the crack path could be approximated by 
connecting a finite number of points on the curve using line segments to create a 
polygonal path (see Figure [T^) . The approximated crack length is computed as 
the sum of the length of the corresponding segments. The actual length could 
be then obtained by computing the limit of the approximated length as the 
segment's length tend to zero. 

When a mesh is introduced to represent the continuum fracture problem 
within the cohesive element formulation, a constraint is introduced into the 
problem due to the inability of the mesh to represent, through its edges (2D) or 
faces (3D), the shape of an arbitrary crack. By following the reasoning in the 
preceding paragraph, a necessary condition for the convergence of the cohesive 
element approach is that the underlying mesh should be able to represent exactly 
a straight line (2D) or plane (3D) oriented in an arbitrary direction as the mesh 
is refined [13]. That is, instead of focusing our attention to the problem of 
representing arbitrary crack shapes, we could concentrate on the ability of the 
mesh to represent a straight line (2D) or a plane (3D), and then think of it as 
a necessary condition to be able to rectificate an arbitrary crack shape. 

In 2-dimensional problems, the ability of a mesh to represent a straight 
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segment is characterized by the path deviation ratio rj, defined as the ratio 
between the shortest path on the mesh edges connecting two nodes, and the 
Euclidian distance between them. The path deviation ratio can be interpreted 
as the maximum relative error in representing a straight line by means of the 
edges of a finite element mesh. Figure [T|d illustrates the relevant quantities to be 
defined for the computation of the path deviation ratio on a given mesh. This 
particular mesh was generated by applying Delaunay's triangulation to a set of 
randomly distributed nodes. The dashed line represents the shortest path over 
the euclidean manifold between the two points marked in bold. The length of 
the dashed line is the euclidean distance between those points, Le- On the same 
figure, the solid line represents the shortest path between the same pair of points 
over the edges of the mesh as computed by means of Dijkstra's algorithm [14] . 
The sum of the lengths of the segments composing the shortest path defines the 
length of shortest path on the graph between those points, Lg. According to 
the previous definitions, the path deviation can be computed as 




(2) 



Other associated measures are the absolute error in representing a straight 
line, given by 77 — 1, or the percentage error given by 100 x (77 — 1). 

Finally, in order to compare path deviation ratios between different kinds of 
mesh and to compare meshes of different sizes, we introduce the non-dimensional 
mesh size, defined as 

N 

A-^f^ (3) 

where N is the total number of edges of the given mesh and hi the length of 
the i*'' edge. That is, the non-dimensional mesh size can be interpreted as the 
inverse of how many average segments are needed to represent a segment of 
length Le- 



2.2 Mesh-induced anisotropy and mesh-induced toughness 

By definition, the path deviation ratio is greater than one, i.e., rj > 1. Con- 
sequently, an arbitrary crack can only be approximated by a discrete curve of 
equal or larger length. For purposes of illustration, let us consider the 4k mesh 
depicted in Figure[2^. The diamond symbols on Figure|3^ show, for such a mesh 
with A w 1/200, how the absolute error changes as function of a chosen direc- 
tion. As observed in the figure, the absolute error (77 — 1) = for the horizontal, 
vertical, and ±45° directions. However, it reaches a peak value near 0.08 for 
intermediate directions. In energetic terms, directions for which the absolute 
error (or path deviation ratio) is lower provide favorable paths for crack prop- 
agation. Clearly, cracks will tend to align along those lower-energy directions 
leading to inaccurate numeric predictions. We term this effect mesh-induced 
anisotropy. 
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Figure 2: (a) typical 4k mesh, and (b) 4k mesh with nodal perturbation. 

In general, when looking at all possible directions at a given point on the 
mesh, the expected value of the path deviation ratio is strictly larger than one, 
i.e., (77) > 1. This, in turn, implies that the expected value of the total energy 
released by the discrete crack {Ef)g is larger than the expected value of the 
energy Ef released by a real crack on the continuum. That is, the introduction 
of a discrete mesh necessarily leads to a larger effective toughness for the discrete 
cohesive model. We call this effect mesh-induced toughness. 

A necessary condition to avoid the undesired mesh-induced anisotropy is for 
the path deviation ratio to be independent on the chosen direction. Additionally, 
a necessary condition to avoid the undesired mesh-induced toughness, is for the 
path deviation ratio to tend to 1 for all possible directions in the mesh as the 
non-dimensional mesh size tend to zero. The satisfaction of these two conditions 
is usually referred to as isoperimetric property. The work of Radin and Sadun 
[15] shows that the pinwheel tiling of the plane has the isoperimetric property. 
Based on this result, Papoulia et al. [13] observed that crack paths obtained 
from pinwheel meshes are more stable as meshes are refined when compared to 
other types of meshes. It is worth noting, however, that pinwheel meshes are 
hard to generate and there is no known extension to the 3-dimensional case. 

Recently, Paulino et al. ^16) introduced two mesh operators that decrease the 
undesired mesh dependent effects observed in 4k meshes. First, they introduced 
a nodal perturbation operator (NP operator) that randomly moves the nodes of 
4k meshes by an amount proportional to the length of the shortest edge in the 
corresponding 4k cell. For example, Figure [^b shows a typical 4k mesh affected 
by the nodal perturbation operator. Subsequently, they introduced a topological 
edge swap operator (ES operator) that allows the originally vertical edges of the 
4k mesh to switch to a horizontal position and vice versa. The idea behind these 
operators is to break the symmetry of the mesh at the geometrical level (NP 
operator) and to relax its structure at the topological level (ES operator) . 
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Figure 3: Error {rj — 1) vs mesh direction for (a) 41c mesii, 4k mesli with nodal 
perturbation, and 4k mesh with nodal perturbation and edge swap, and (b) 4k mesh 
with nodal perturbation and edge swap, and random mesh. 



As noted in their work, both operators produce an increment in the value of 
the path deviation ratio in the directions that showed no error for the 4k mesh, 
while at the same time they reduce it for the intermediate directions. They 
also show that the reduction of the path deviation ratio in the intermediate 
directions is larger than the increment produced in the the horizontal, vertical, 
and ±45° directions. The combination of these effects results in an effective 
reduction of the mean value of the path deviation ratio (ry) . It is worth noting, 
however, that even though the mesh induced-toughness is reduced in this way 
(that is, the mean value of the path deviation ratio is reduced), the meshes 
under consideration still exhibit high anisotropy, only evident when considering 
the full polar plot of the absolute error (or path deviation ratio). We illustrate 
this effect by analyzing a 4k mesh with nodal perturbation factor NP/ = 0.3 (4k 
with NP in Figure |3K) , and a 4k mesh perturbed by the same amount with the 
extra possibility of swapping edges (4k with NP+ES in Figure [5^). As Figure 
[3^ shows, for meshes with non-dimensional mesh size A « 1/200, the maximum 
attained value of the path deviation ratio is decreased by the application of the 
NP operator. Furthermore, further application of the ES operator results in 
even lower values of the path deviation ratio. However, for all three cases the 
path deviation ratio highly depends on the direction under consideration. That 
is, both the NP and ES operators seem to have little effect on the mesh-induced 
anisotropy. 

To conclude this section, we observe an interesting feature that arises from 
analyzing the polar behavior of a random mesh, i.e., a mesh generated through 
the application of Delaunay's triangulation algorithm to a set of randomly placed 
points. Figure depicts the absolute error dependence on the direction for a 
random mesh and a 4K mesh subject to the NP and ES operators. The adopted 
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value of the nodal perturbation factor is NP/ — 0.3, and both meshes have 
comparable values of A « 1/200. It is interesting to note that whereas the 4k 
mesh exhibits a lower mean value of the absolute error ratio {{rj — 1) « 0.37 
vs (?7 — 1) ~ 0.43), the random mesh shows a much better isotropic behavior 
(std(?7 — 1) sa 0.002 vs std(?7 — 1) w 0.016). Needles to say, ideally we would like 
to have a mesh that reflects a compromise between these two behaviors, that is, 
a low absolute error that remains constant for all possible directions. 

3 K-means and conjugate-directions meshes 
3.1 Isotropy through K-means meshes 

Based on the previous observations, we wish to generate random meshes in 
order to achieve isotropy in the sense of the path deviation ratio. At first, 
purely random meshes look attractive: even though they generate path deviation 
ratios with mean value slightly higher than those obtained by 4k meshes under 
the effects of NP and ES operators, they have the advantage of not exhibiting 
any undesired anisotropy. However, purely random meshes are very irregular 
and present highly distorted triangles rendering them unfit for finite element 
computations. To address this issue, we propose to smooth purely random 
meshes by subjecting them to a K-means clustering algorithm. The K-means 
clustering algorithm partitions spn ■ n observations into n clusters in which each 
observation belongs to the cluster with the nearest mean. Then, a so called 
K-means mesh can be obtained as follows: 

• spn ■ n uniformly random points are clustered into n nodes. 

• The resulting nodes are then triangulated by means of Delaunay's algo- 
rithm. 

For the purposes of this work, K-means clustering was implemented through 
Lloyd's algorithm. For more details on Lloyd's algorithm and the relationship 
between K-means clustering and centroidal Voronoi tesselations, see [17]. Figure 
m shows a typical K-means mesh obtained obtained by clustering 128 • 512 = 
65536 random points into 512 nodes {spn — 128, n — 512). 

To better illustrate the effect of the K-means clustering algorithm on random 
meshes, we demonstrate the infiuence of the clustering parameter spn on both 
the spatial node distribution, and on the triangle quality distribution of the 
resulting meshes. 

Figure [5] shows the effect of the clustering parameter spn on the spatial 
distribution of mesh nodes and the corresponding Fourier spectrum. The Figure 
is organized in a matrix layout where rows (from top to bottom) correspond to 
nodes clustered with spn values of 1, 8, 32 and 128. The left column shows the 
spatial distribution of nodes, the center column the corresponding 2D Fourier 
power spectrum, and the right column the radially averaged power spectrum. 
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Figure 4: Typical K-means mesh obtained by clustering 128 ■ 512 = 65536 random 
points into 512 nodes {spn = 128, n = 512). 

The spn = 1 case corresponds to a uniformly random distribution of nodes. 
The corresponding 2D- Fourier power spectrum exhibits an even frequency con- 
tent typical of white noise. As the spn parameter increases, nodes tend to 
distribute more evenly in space, and the corresponding 2D-Fourier spectrum in- 
dicates that meshes remain isotropic (no radial variation in the spectrum) after 
the clustering algorithm. Due to this angular invariance, all spectral informa- 
tion can be exhamined through the radially averaged power spectrum. As it can 
be observed from the right column on Figure [SJ the clustering algorithm filters 
the lower frequencies in the Fourier space, resulting in a typical blue noise power 
spectrum. That is, the lower the value of spn the closer the tendency to generate 
to a purely white noise mesh, where as spn increases the mesh approaches the 
behavior of a blue-noise mesh, which preserves isotropy and randomness while 
filtering low frequencies resulting in more even spacing between nodes. 

Figure ini shows the effect of the clustering parameter spn on the mesh quality 
distribution. The Figure is organized in a matrix layout where rows (from top to 
bottom) correspond to nodes clustered with spn values of 1, 8, 32 and 128. The 
left column shows the obtained mesh, and the right column the distribution of 
the mesh quality parameter q — 2r/R within the mesh, where r is the inradius 
and R is the circumradius of a triangle in the mesh. The chosen mesh quality 
parameter can take values between and 1 , indicating a sliver (or zero- volume 
triangle) and 1 indicating an equilateral triangle. Clearly, as the value of spn 
increases the corresponding mesh quality distribution improves drastically. As 
a result, K-means meshes with higher value of the cluster parameter spn tend 
to contain a larger proportion of nearly equilateral triangles. 

To conclude our study of K-means meshes, we analyze the effect of the 
clustering parameter spn on the metrics relevant to our original problem: the 
path deviation ratio. Figure [7^ depicts the value of the mean value of the 
path deviation ratio as a function of spn for a mesh with non-dimensional mesh 
size A « 1/400. As it can be observed, an increment on the value of spn 
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Figure 5: Influence of the smoothing parameter spn on mesh node distribution. Spa- 
tial node distribution (left column), 2D-Fourier power spectrum (center column), and 
radiaUy averaged power spectrum (right column). 
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Figure 6: Influence of the smoothing parameter spn on mesh quality. Mesh (left 
column) and corresponding mesh quality distribution (right column). 
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Figure 7: Influence of the smoothing parameter spn on (a) mean(r;), and (b) std(77) 
for a K-means mesh. 



initially produces a slight reduction on the mean value of the path deviation 
ratio. After spn is further increased, the mean value of the path deviation ratio 
slightly increases, saturating at a value similar to that observed for the spn = 1 
case. Additionally, Figure [TId indicates that the standard deviation of the path 
deviation ratio is not much affected by the value of spn. That is, the smoothing 
process performed through the clustering algorithm improves the mesh quality 
while not affecting the behavior of path deviation ratio. 

In summary, K-means meshes provide a good alternative to 4k meshes with 
NP and ES operators, as they provide slightly higher values for the mean value 
of the path deviation ratio while exhibiting a perfectly isotropic behavior. 



3.2 Reducing mesh induced toughness through conjugate- 
directions meshes 

The good isotropic behavior of the K-means meshes could be attributed to the 
ability of the mesh, at any node, to provide a set of random directions for the 
crack to continue its propagation towards a neighboring node. Similarly, the 
mean value of the path deviation ratio can be reduced by enriching the space of 
possible directions available at a point for a crack to propagate to the next point 
in the mesh, as done by the ES operator in 4k meshes. Towards this end, we 
introduce a new type of mesh, which we term conjugate-directions mesh. This 
new family of meshes preserves isotropy by using as starting point an existing 
K-means mesh, and enriches the space of possible directions by introducing 
topological modifications to it. Conjugate-directions meshes are obtained as 
follows: 

• A K-means mesh is generated 
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• A barycentric subdivision is applied to the initial K-means mesh 



An algorithm for performing a barycentric subdivision to simplicial meshes 
of any dimension is detailed in , and explained in the subsequent paragraphs 
for completeness. Towards this end, we first describe (briefly) the concept of 
simplicial complex and subdivisions of a simplicial complex (triangulation) . 

A cell complex is a collection of objects, or cells, to which a precise dimen- 
sion can be assigned. At its fundamental level, a 3-dimensional object can be 
represented by a collection of 0-ceUs (vertices), 1-cells (edges), 2-cells (faces) 
and 3-cells (volumes). In the subsequent paragraphs, p-dimensional cells are 
denoted as Cp. In addition, Ep{S) represents the collection of p-dimensional 
cells in a complex S. The dimension dim 5 of the cell complex S is the largest 
dimension of any of its cells. Special classes of cell complexes are obtained when 
the cells are restricted to be of a certain type. 

If all cells in a cell complex are simplices, the cell complex is said to be 
a simplicial complex. Simplicial complexes arise naturally as a result of the 
triangulation of solids. Specifically, a simplicial complex S in M" is a collection 
of simplices in R" such that (i) every face of a simplex of S is in S, and (ii) the 
intersection of any two simplices of 5 is a face of each of them. 

We recall that the simplex a spanned by a geometrically independent point 
set vq, Vp in M" is 



a=lxeR": = 1, ^ A^w^ = x, A^ > Vz = 0, 1, ^ (4) 



where the numbers Xi are the barycentric coordinates of x with respect to 
Vq, Vp. The barycentric coordinates of a point are uniquely determined by it. 
The points fo, Vp are the vertices of the simplex a and p its dimension. Any 
simplex generated by a subset of vq, ...,Vp is a proper face of a and the union of 
all proper faces is the boundary of the simplex. The notation -< ec, signifies 
that is a face of Cq. 

Let 5 be a cell complex. A subdivision complex of S is obtained by sub- 
dividing its cells into finer cells. More precisely, a complex S* is said to be a 
subdivision of S if (i) every cell of S* is contained in a cell of S, and (ii) every 
cell of S is the union of finitely many cells of S* . These conditions particularly 
imply that the union of the cells of S* equals the union of the cells of S and, 
hence, iiS*| — \S\ as sets. The way in which a subdivision complex S* is nested 
within the supercomplex S may be described by means of an inclusion map 
f : S* ^ S. This map assigns to every cell ep of S* the cell of S that 
contains it. 

Barycentric subdivision generates a particular class of subdivision complexes. 
Thus, the barycentric subdivision of a simplicial complex is obtained by means 
of a uniform refinement. If Up — [wo, Vp] is a p-simplex of a simplicial complex 
S, its barycenter is the point 



p 



p 




(5) 



12 



Figure 8: Barycentric subdivision of a triangle. 



i. e., the barycenter is the point in the interior of that has equal barycentric 
coordinates. Given a simplicial complex S its barycentric subdivision S* consists 
of all simplices [ao, dp] such that ctq >- ... >- d'p. The details of the baycentric 
subdivision algorithm are depicted in Algorithm [1] for completeness. As an 
example, Figure [8] shows the barycentric subdivision of a triangle. It is worth 
noting that the barycentric subdivision of a triangle will produce 6 new triangles 
whose aspect ratio is worse than that of the the original one. 



Algorithm 1: Barycentric subdivision. 



Input: simplicial complex S, n = dimS. 
Output: barycentric complex S* . 

/* Initialize barycentric complex S* and inclusion map / */ 
1 for {p — 0, n) do 

2 
3 
4 



for (cp e Ep{S)) do 
insert v in Eq{S*); 
insert {cp, v) in /; 



/* Compute connectivity table C'T{S) */ 

5 for (e„ e En{S)) do 

6 1^ insert [e„] in CT{S); 

7 for {p = n, 1) do 



8 
9 
10 



for ([e„,...,ep] G CT{S)) do 
for (cp-i -< Bp) do 
|_ append Cp-i to [e„,...,ep]; 



/* Compute connectivity table CT(S*) */ 

11 for ([e„,...,eo] G CT{S)) do 

12 \_ insert [/(e„), /(eo)] in 

13 construct the simplicial complex S* from CT{S*); 

14 return S* 



An interesting feature of the barycentric subdivision, and the reason why we 
chose it to generate conjugate-directions meshes, is that in the limiting case of 
equilateral triangles, the barycentric subdivision of a triangle provides an orthog- 
onal direction for each direction provided by the edges of the original triangula- 
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(a) 



Figure 9: Typical conjugate-directions mesii obtained from a K-means mesli witli 
spn = 128 and n — 512. 

tion. This can be clearly observed in Figure [51 where there is a newly introduced 
dashed line perpendicular to each edge of the original triangle. In this way, for 
K-means meshes with high value of spn (for which we already illustrated their 
higher content of equilateral triangles), their barycentric subdivision enriches 
them with a set of directions that are almost orthogonal to the existing edges. 
Figure [9] shows a typical conjugate-directions mesh obtained from a K-means 
mesh with spn — 128 and n — 512. 

Figure [10] shows the values of the absolute error vs direction for a 4k mesh 
with NP and ES, a K-means mesh, and a conjugate-directions mesh. All 
meshes are comparable in terms of their non-dimensional size A « 1/200. The 
conjugate-directions mesh is not only isotropic, but the mean value of the ab- 
solute error is significantly lower 0.018) than that observed for the other 
meshes (« 0.04). 

It is worth noting that, since conjugate-directions meshes are derived from 
K-means meshes, the spn parameter of the seed mesh will play a role in the 
performance of the newly generated one. Figure [TT] shows the effect of the spn 
parameter on the path deviation ratio for a conjugate-directions mesh with A w 
1/400. As shown in lllb . an increment of spn produces a consistent reduction 
on the mean value of the path deviation ratio. For spn = 512 the mean value 
of the path deviation ratio reaches a value (77) « 1.015, approximately half 
(in error terms) of that observed for the barycentric mesh obtained from a 
random {spn — 1) case. Additionally, Figure [TTb indicates that the standard 
deviation of the path deviation ratio decreases as spn increases. This implies 
that, when a K-means mesh with high value of spn is adopted as a generator of a 
conjugate-directions mesh, the resulting triangulation produces a more compact 
distribution of the path deviation ratio. 
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Figure 10: Error {-q — 1) vs mesh direction for a 4k mesh with NP and ES, a K-means 
mesh, and a conjugate-directions mesh. All non-dimensional mesh sizes correspond to 
A ^ 1/200. 

In summary, conjugate-directions meshes provide an even better alternative 
to 4k meshes with NP and ES, as they provide significantly better values for 
the mean value of the path deviation ratio while being perfectly isotropic. It 
is worth noting that the original seed mesh must be of good quality since the 
barycentric subdivision stage in the generation of a conjugate-directions mesh 
produces a deterioration in the quality of the triangulation. 

3.3 Convergence of K-means and conjugate-direction meshes 

To conclude this section, we perform a numerical convergence analysis of 4k, 
K-means and conjugate-direction meshes for the non-dimensional mesh size A 
tending to zero. All K-means and conjugate-directions meshes considered in 
this numerical study of convergence were generated with a clustering parameter 
spn = 128. The nodal perturbation factor adopted for all 4k meshes is NP/ = 
0.3. 

Results show that convergence of K-means meshes, in the sense of the mean 
value of the path deviation ratio, is similar to that of 4k meshes with NP as the 
non-dimensional mesh size is reduced, see Figure 112b .. Overall, 4k meshes with 
NP and ES exhibit lower mean values of the path deviation ratio within the 
mesh size range under consideration. This difference seems to decrease as the 
meshes are refined. Numerical evidence also shows that the mean value of the 
path deviation ratio tends to saturate around 1.04 for 4k with NP and K-means 
meshes, and around 1.036 for 4k mesh with NP and ES. That is, these three 
types of mesh have an intrinsic roughness that prevents them from representing 
straight lines with errors below 4% or 3.6% respectively no matter how refined 
they are. 

On the other hand, the standard deviation of the path deviation ratio for 
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Figure 11: Influence of the smoothing parameter spn on (a) mean(77), and (b) std(r;) 
for a conjugate-directions mesh. 



K-means meshes decreases at a fastest rate when compared to both types of 4k 
meshes, see Figure [T2b . This indicates that even though both kinds of mesh 
exhibit similar mean values of the path deviation ratio for all mesh sizes, K- 
means meshes show less dispersion as the mesh is refined, while 4k meshes 
saturate very quickly. The saturation observed for both types of 4k mesh is an 
indicator that, in these cases, the mesh-induced anisotropy does not disappear 
as the meshes are refined. 

In addition, the same numerical experiment shows no indication of saturation 
for conjugate-directions meshes in the studied range of non-dimensional mesh 
sizes. Furthermore, the mean value of the path deviation ratio is significantly 
smaller (on the order of 1/3 to 1/2 in absolute error terms) for conjugate- 
directions meshes when compared to the other meshes for all non-dimensional 
mesh sizes considered in this study. The standard deviation of the path deviation 
ratio seems to converge at a slower rate for this kind of mesh when compared 
to K-means meshes, but at a faster rate when compared to 4k meshes. It is 
worth noting, however, that it seems that the convergence rate on the standard 
deviation of the path deviation ratio tends to accelerate as the mesh is refined. 



4 Summary and final remarks 

We illustrated how the introduction of a discretization produces mesh-induced 
anisotropy and mesh-induced toughness when dealing with cohesive element 
modeling of arbitrary cracks. Subsequently, we analyzed those effects through 
the introduction of polar plots of the path deviation ratio, and showed that 
randomness plays a crucial role in the reduction of the former effect. 

We observed that 4k meshes with NP and ES provide a better alternative 
to plain 4k meshes in the sense that they reduce the mesh-induced toughness. 
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Figure 12: Numerical convergence analysis of 4k, K-means, and congugate-direction 
meshes. Mean value of path deviation ratio (a) and standard deviation of path devi- 
ation ratio (b) as a function of the non-dimensional mesh size. 



However, when looking at the full polar plot of the path deviation ratio vs 
mesh direction, it becomes evident that mesh-induced anisotropy, while reduced, 
still persists. That is, regardless of the type of 4k mesh (with or without NP 
and/or ES), there is always an associated anisotropy to it. In addition, this 
effect seemingly does not disappear as the mesh is refined. This has profound 
implications on the interpretation of mesh refinement algorithms involving 4k 
meshes: at some point, there is no gain from mesh refinement no matter how 
much these meshes are refined. 

On the other hand, K-means meshes were shown to be a good compromise 
between isotropy and good element quality. They are isotropic thus not provid- 
ing preferred crack propagation directions for any mesh size. However, K-means 
meshes saturate, in the sense of the mean value of the path deviation ratio, as 
they are refined. This saturation results in the introduction of mesh-induced 
toughness even in the limiting case of non-dimensional mesh size tending to 
zero. A good feature of K-means meshes is that they do not seem to saturate 
in the sense of the standard deviation of the path deviation ratio, which implies 
that the error concentrates around its mean as the mesh is refined. This would 
allow to adopt K-means meshes in cohesive element modeling as long as the 
proper correction term is introduced in the toughness of the material. 

Finally, we proposed a new type of mesh, termed conjugate-directions mesh, 
which significantly reduces the undesired mesh-dependent effects for meshes 
of practical size. Our approach combines the use of barycentric subdivision 
with K-means meshes to provide a greater distribution of directions for fail- 
ure surfaces to propagate. Numerical evidence suggest that, due to the exhib- 
ited isotropy and considerably lower mesh-induced toughness effect, conjugate- 
direction meshes could be good candidates for cohesive element analysis of crack 
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propagation problems where crack paths are not known a priori. Our hope is 
that, by designing meshes to mitigate the effects of mesh topology, we can invoke 
adaptive insertion along pre-defined element boundaries to achieve arbitrary 
crack propagation. The great promise of this research effort is that arbitrary 
crack propagation could be achieved through mesh design (pre-processing) . 
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